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We compute rates and pathways for nucleation in a sheared two dimensional Ising model with 
Metropolis spin flip dynamics, using Forward Flux Sampling (FFS). We flnd a peak in the nucleation 
rate at intermediate shear rate. We analyse the origin of this peak using modified shear algorithms 
and committor analysis. We find that the peak arises from an interplay between three shear-mediated 
effects: shear-enhanced cluster growth, cluster coalescence and cluster breakup. Our results show 
that complex nucleation behaviour can be found even in a simple driven model system. This work 
also demonstrates the use of FFS for simulating rare events, including nucleation, in nonequilibrium 
systems. 



INTRODUCTION 



The nucleation of a stable phase from a metastable 
one is a ubiquitous and important phenomenon. Most 
progress in understanding the physics of nucleation has 
been made for "quasi-equilibrium" systems, in which the 
system dynamics obeys detailed balance and the transi- 
tion is from a metastable to a thermodynamically stable 
state. However, many important nucleation processes 
both in nature and in industry happen in driven sys- 
tems, such as those under shear, whose dynamics do not 
obey detailed balance. Despite its importance, nucle- 
ation in driven systems remains poorly understood. In 
this paper, we compute rates and transition paths for a 
driven nucleation process: nucleation under shear in a 
two dimensional Ising model. We use the recently de- 
veloped Forward Flux Sampling rare event simulation 
method P, We find that the nucleation rate shows 

a striking nonmonotonic dependence on the shear rate, 
and that this is due to an interplay between three shear- 
mediated effects: shear-enhanced cluster growth, cluster 
coalescence and cluster breakup. 

Nucleation under shear remains poorly understood 
[1, 0|. It is expected that high enough shear rates will 
impede nucleation. Some studies of crystal nucleation 
0, 0, S] find that nucleation rates decrease monotoni- 
cally with shear rate; others suggest that crystallisation 
may be enhanced at low shear rates [l^, [ill, [l2|- A 
recent experimental study found a minimum in the crys- 
tal nucleation rate as a function of shear rate for charged 
colloids [l^. Crystallisation from sheared glassy states is 
even more complicated, both experimentally and numer- 
ically [l3, [l3|. For binary mixtures [lB| and isotropic- 
to-lamellar transitions shear is predicted to increase 
the critical temperature. Physical mechanisms for the 
effect of shear on nucleation may include hydrodynamic 
effects, cluster coalescence, cluster breakup, layering due 
to the shear fiow, and suppression of polydispersity. In 



this work, we study an idealised model in which many of 
these effects are not included (perhaps most significantly, 
transport processes are not modelled). Our motivation is 
to provide a fundamental basis on which to build an un- 
derstanding of more complex systems. Our results may 
however be relevant to driven solid materials (itI . [18| . 

The Ising model provides a paradigm for many phe- 
nomena in statistical physics, including nucleation. Nu- 
cleation in this model, in the absence of external driving 
has been extensively studied [H, [13, [U, [13, [13, H El, 0, 
[23,[13|- Ising models have proved a valuable tool for test- 
ing the Classical Nucleation Theory (CNT)[2i,|33,|3il, in 
which nucleation is coarse-grained to one dimension, the 
largest cluster size coordinate, and modelled as a transi- 
tion over a free energy barrier that arises due to compe- 
tition between the favourable chemical potential of the 
growing cluster and its unfavourable interfacial free en- 
ergy. An attempt has been made to extend the CNT 
to sheared systems [111. In the absence of shear, tran- 
sition path analysis has shown the importance of order 
parameters other than the largest cluster size in the nu- 
cleation mechanism, in both two and three dimensions 
[33, [131 • This large body of information on nucleation 
in the absence of driving makes the Ising model an at- 
tractive test system for nonequilibrium nucleation prob- 
lems. Metastability and nucleation of nonequilibrium 
steady states generated by coupling to two different heat 
baths has been studied in a two-dimensional Ising model 
[H, [111- Although, to our knowledge, nucleation under 
shear has not been investigated for the Ising model, the 
maximum likelihood path has recently been computed for 
nucleation under shear in a finite system defined by a dy- 
namical equation for the nucleation order parameter, in 
the absence of applied field |33|. In this paper, we study a 
sheared two-dimensional Ising model. We find that even 
for this highly simplified system, nucleation under shear 
is a complex process that depends on multiple physical 
mechanisms. 

The Forward Flux Sampling (FFS) method used in this 
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work allows the computation of rate constants, transition 
paths and stationary probability distributions for rare 
events in equilibrium or nonequilibrium systems 0,0, [s^]. 
Rare events, such as nucleation, are notoriously difficult 
to simulate, because the waiting time between events is 
typically much longer than the timescale of the event 
itself, meaning that few, if any, events are observed in 
a typical simulation run. Rare event simulation meth- 
ods developed for equilibrium systems include Bennett- 
Chandler methods [3§|, transition path sampling [4o| . 
(partial path) transition interface sampHng [4l|, mile- 
stoning and string methods [Hi. Both TPS and 
the string method have been applied to Ising nucleation 
[s^, 113, Ell . However, these methods require knowledge 
of the steady state phase space density, making them un- 
suitable for nonequilibrium problems. FFS does not re- 
quire knowledge of the phase space density. The method 
uses a series of interfaces in phase space between the ini- 
tial and final states, defined by an order parameter which 
need not be the reaction coordinate. In earlier work, we 
and others have shown that FFS correctly reproduces the 
nucleation behaviour of a two dimensional Ising model in 
the absence of shear [2^ . [38| . 

In section [ill we give details of our simulation model 
and the FFS method applied to this system. In section 
lllli we present results for the nucleation rate as a function 
of the shear rate. We then analyse the physical mecha- 
nism behind the suppression of nucleation at high shear 
rates in section Uvl In section |Vl we discuss the roles 
of shear-enhanced cluster growth and coalescence in the 
enhancement of nucleation at low shear rates. We test 
our ideas with an analysis of cluster growth in section 
IVll with a comparison to a modified shear algorithm in 
section IVIIl and with an analysis of the transition state 
ensemble in section IVIIII Finally, we present our conclu- 
sions in section HXl 

II. SIMULATION DETAILS 

The sheared two-dimensional Ising model 

Our system consists of a two-dimensional Lx L square 
lattice of up-down spins with nearest-neighbour spin-spin 
interactions, coupling to an external magnetic field, and 
periodic boundary conditions in the x and y directions. 
The Hamiltonian for the spin-spin and spin-field interac- 
tions is 

H = -j'^aiffj - h^^ai, (1) 

ij i 

where ai — ±1 is the state of spin i, J is the coupHng 
constant between neighbouring spins and h the external 
magnetic field - both in units of fc^T. The prime in- 
dicates that the sum is restricted to nearest-neighbour 
interactions. Our simulations use Metropolis spin-fiip 
dynamics. In each Monte Carlo cycle, we make L x L 



attempts to fiip a spin. In each attempt, we choose a 
spin at random, attempt to fiip it, and accept or reject 
the fiip according to the Metropolis rule. An alterna- 
tive choice of dynamics, not considered here, would be 
the Kawasaki scheme in which up spins diffuse between 
lattice sites f4$|. 

All the simulations described here use a lattice of size 
L = 65, and coupling constant J = 0.65fc_BT. We apply 
an external magnetic field h = O.OSfcsT in all simula- 
tions. Our coupling constant J is larger than the critical 
coupling Jc « 0.44fcsT |4fl]. Considering the system in 
the absence of shear, the thermodynamically stable state 
is ferromagnetic, with net positive magnetisation, mean- 
ing that the system tends to have the majority of its spins 
in the up state (cr = +1). The alternative ferromagnetic 
state with net negative magnetisation (most spins in the 
down state) is metastable, and if initiated in a predom- 
inantly down state, the system will remain in that state 
for a significant time before undergoing a nucleation tran- 
sition to the thermodynamically stable up state [1^. We 
are interested in the rates and pathways for this transi- 
tion. In the absence of shear, this system is identical to 
that investigated by Sear [Hi and previously by some of 
us [Hi, except that we now use a larger box size, since we 
have found that the nucleation rate is more sensitive to 
system size in the presence of shear. For the shear rates 
used in this study, L = 65 is large enough to ensure that 
our computed nucleation rate is independent of L. The 
free energy barrier in the absence of shear is ~ 22k bT ^ 
so that we are working at moderate supersaturation. We 
therefore expect nucleation to proceed via the growth of 
a single large cluster of up spins. 

We apply shear to the system using a method similar 
to that of Cirillo et al After each Monte Carlo cy- 
cle, we make x L attempts to shear the system {Ms 
is the number of attempted shear steps per row per MC 
cycle). In each attempt, we carry out a shear step with 
probability Ps- A shear step consists of choosing a row 
js at random, and shifting all lattice sites with j > js to 
the right by one lattice site. The net result is that row j 
is shifted to the right by on average jMgPs lattice sites 
per Monte Carlo cycle. The shear rate is thus given by 
7 = MgPs- We have verified that the choice of Mg and 
Ps does not matter for a given product 7. We note that 
care must be taken to maintain the correct identity of 
the neighbour sites in the periodic image lattices above 
and below the simulation box - after a shifting move, 
the identity of these neighbours is changed. Our method 
for achieving this is described in Appendix [Aj This al- 
gorithm imposes, on average, a linear velocity field on 
the underlying lattice. In a real physical system, the ve- 
locity field is not imposed externally but emerges as a 
consequence of the internal dynamics of the system, so 
that the shear algorithm used here is somewhat artificial. 
However, our purpose here is to investigate the effects of 
a linear velocity field on the system; furthermore, this al- 
gorithm has the advantages of being simple to implement 
and homogeneous across the simulation box. 
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Forward Flux Sampling 

We have used the Forward Flux SampHng (FFS) 
method [J, [3, S| to calculate nucleation rates and tran- 
sition paths for the formation of the steady state with 
predominantly up spins (the up state), from the steady 
state with predominantly down spins (the down state). 
This rare event sampHng method uses a series of inter- 
faces in phase space between the initial and final states 
to force the system from the initial state A to the final 
state -B in a ratchet-like manner. An order parameter 
X{x) is defined (where x represents the phase space co- 
ordinates), such that the system is in state A if A(a;) < Aq, 
and it is in state B if A(x) > A„, while a series of non- 
intersecting planes (interfaces) Xi {0 < i < n) lie be- 
tween states A and B, such that any path from A to B 
must cross each interface, without reaching A^+i before 
Xi. Provided enough configurations are obtained at the 
first interface to ensure good sampHng, the choice of or- 
der parameter A(a;) should not affect the calculated rate 
constant or transition paths {i.e. X{x) need not be the 
true reaction co-ordinate) - although it may affect the 
computational efficiency of the method. 

Full details of the FFS method are given in Refs 0,[3|, 
and a detailed analysis of its computational efficiency is 
given in Ref Briefly, the transition rate / from A to 
B is decomposed as [4ll.l48l|: 

n-l 

I = $A,n = ^A,oP{Xn\Xo) = $A,0 J] Pi^^+l\>^^)■ (2) 

where $A,n is the average flux of trajectories crossing 
from A to B, ^a.o is the average flux of trajectories cross- 
ing Ao in the direction of B, P(A„|Ao) is the probability 
that a trajectory that crosses Aq in the direction of B will 
eventually reach B before returning to A, and P{Xi+i\Xi) 
is the probability that a trajectory which reaches Xi, 
having come from A, will reach Ai+i before returning 
to A. The flux $a,o is computed using a simulation in 
the A state, during which conflgurations corresponding 
to crossings of the flrst interface Ao coming from A are 
also stored. This collection of conflgurations is then used 
to initiate trial runs which either reach the next interface 
Ai or go back to Ao, generating an estimate for the con- 
ditional probability P(Ai|Ao) as well as a new collection 
of conflgurations at Ai; the trial run procedure is then 
iterated until B is reached. The rate constant / is then 
obtained from Eq.(l2]), and a correctly weighted collection 
of trajectories from A to B is obtained by tracing trial 
runs that eventually arrive at A„ , via successive interfaces 
back to A. In practice, rather than storing all conflgu- 
rations for all trial runs during the FFS procedure, it is 
sufficient to store the initial collection of configurations 
at Ao , together with limited information about each con- 
figuration in each collection at intermediate interfaces Ai, 
indicating its "parent" configuration in the collection at 
the previous interface Ai_i, and the value of the random 
number seed used to initiate the relevant trial run. In 



this way, transition paths can easily be reconstructed af- 
ter the FFS sampling procedure, without the need for ex- 
cessive data storage. In an earlier study, we have shown 
that FFS correctly reproduces the nucleation behaviour 
of a two dimensional Ising model in the absence of shear 

[ai- 

For the simulations described in this paper, the param- 
eter A was defined as the total number of up spins in the 
simulation box. This is a global order parameter: an al- 
ternative might be to use the size of the largest cluster of 
up spins. For FFS, we do not expect the rate constant or 
the transition path ensemble to depend on the choice of 
order parameter (we will later analyse trajectories using 
the largest cluster size). Others have experienced sam- 
pling problems when using global order parameters in 
FFS [43|; this was not the case in this work. We used 39 
interfaces for our FFS calculations (except in section [Vl, 
defining the A state at A < Aq where Ao = 25 up spins, 
and the B state at A > A„ where A„ — 2005 up spins (the 
total number of spins in our system being 65x65 = 4225) . 
The spacing between interfaces varies between 5 and 200 
up spins. We collected 1000 configurations at interface Aq 
and repeated each FFS calculation 25 times. The num- 
ber of trials at each interface varied between 1300 and 
7000. Our results do not depend on the precise choice of 
the number or position of the interfaces. 



III. NUCLEATION RATE AS A FUNCTION OF 
SHEAR RATE 




Y 

Figure 1: Rate of homogeneous nucleation 7 as a function of 
shear rate 7 for ft = O.OSfcsT and J = 0.65fcsT. 



Figure [T] shows the rate / of homogeneous nucleation 
as a function of the shear rate 7. The nucleation rate 
shows a striking nonmonotonic dependence on 7. For 
low shear rates, / increases apparently linearly with 7, 
before reaching a maximum around 7 = 0.06. For shear 
rates 7 > 0.06, / decreases nonhnearly with increasing 
7. For 7 = 0, our result is in good agreement with the 



4 



value of 3.3 x 10 per MC cycle per site obtained by 
Sear @. 
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Figure 2: Configurations from the transition state ensemble, 
obtained as described in Appendix [B] (a): 7 = 0.0, (b): 7 — 
0.06, (c): 7 = 0.12. 



Figure [2] shows representatives of the transition state 
ensemble (TSE; these are transition path configurations 
from which a newly initiated trajectory has probability 
Pb = 0.5 of reaching B before A) for shear rates 7 = 0.0, 
7 = 0.06 and 7 = 0.12. It is clear that the shape of 
the growing cluster is strongly affected by the shear. An 
observation of the transition paths shows that, for high 
shear rates, the growing cluster eventually connects with 
its periodic images to form a horizontal stripe across the 
box, which then expands vertically to fill the box. This 
has also been seen for nucleation under shear in the un- 
physical case of no supersaturation [Hi. In our simula- 
tions, stripe formation occurs well beyond the transition 
state, and does not affect the nucleation rate (since we 
have verified that / is independent of L) . 

In the following sections, we attempt to elucidate the 
physical origin of the nonmonotonic dependence of / on 
7 shown in Figure [H We first consider the origin of the 
decrease in / with 7 at high 7, and then turn to the 
mechanisms behind the increase in 1(7) for low 7. 



IV. SUPPRESSION OF NUCLEATION AT HIGH 
SHEAR RATE 

We first seek an explanation for the decrease in nu- 
cleation rate / with shear rate 7 for 7 > 0.06 in Fig- 
ure [H Figure [2] shows that the growing clusters become 
elongated in the direction of the shear. The extent of 
this elongation is governed by a balance between the fre- 
quency of shear steps and the growth rate of the cluster. 
It seems intuitive that for high shear rates, the elonga- 
tion due to the shear will exceed the rate at which the 
cluster can grow, leading to shear-induced breakup of the 



growing cluster, and a corresponding decrease in the ho- 
mogeneous nucleation rate. 
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Figure 3: Configurations from the transition state ensemble 
with rattle shear, obtained as described in Appendix [B] (a): 
7 = 0.0, (b): 7 = 0.06, (c): 7 = 0.12. 




Figure 4: / versus 7 for h — 0.05kBT, for regular shear (cir- 
cles) and rattle shear (squares). The regular shear results are 
the same as in Figure [T] 

To test this hypothesis, we performed a set of simu- 
lations in which the direction of the shear (to the right 
or to the left along the x axis) was chosen at random 
for each shear step. We call this algorithm "rattle shear". 
On average, the system makes as many row shifts to the 
right as it does to the left, so we do not expect clusters 
to be elongated by the shear. This is confirmed in Figure 
[31 which shows that TSE configurations are not notice- 
ably elongated, even for high shear rates. Fig [4] shows 
/ versus 7 for the rattle shear (squares), as well as for 
the "regular" shear algorithm (circles). As expected, the 
regime in which / decreases with 7 has been abolished 
for the rattle shear algorithm, at least within this range 
of 7 values. This appears to confirm our hypothesis that 
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the decrease in / for large 7 is due to shear-induced elon- 
gation of the growing cluster, leading to eventual cluster 
breakup. 

Figure |4] also shows that nucleation is enhanced less 
strongly at low shear rates for rattle shear than for the 
regular shear algorithm. This suggests that multiple 
shear steps occurring in the same direction are an impor- 
tant feature in the nucleation enhancement mechanism; 
a topic which is explored further in the next sections. 

V. ENHANCEMENT OF NUCLEATION AT 
LOW SHEAR RATE 
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Figure 5: A square cluster (left) undergoes one shear step 
to make a shape with two concave kinks and two additional 
corner sites. Kink sites are labelled K and corner sites C. 

We turn next to the increase in nucleation rate / with 
shear rate 7 observed in Figure [T] for low shear rates, 
7 < 0.06. This behaviour could have (at least) two pos- 
sible origins. Firstly, the shear algorithm changes the 
shape of the growing cluster and this may increase its 
tendency to grow. One way in which this could happen 
is by increasing the surface roughness of the cluster: our 
shear algorithm creates kinks in the growing cluster, and 
these are favourable sites for further cluster growth. The 
rate of growth due to spin flips is therefore likely to be 
enhanced by the shear. This mechanism is illustrated in 
Figure HJ where two kinks (labelled K) and two corner 
sites (labelled C) are created by a shear step. Although 
the corner sites have a tendency to flip to the down state, 
this is counteracted by the greater tendency (due to the 
applied field) of the kink sites to fiip to the up state. 




Figure 6: The shear algorithm can drive isolated spins to- 
gether into small clusters. The top three rows are shifted to 
the right by a shear step, causing two isolated up spins to fuse 
into a small cluster. 



(b) 



(d) 




Figure 7: The shear algorithm can drive isolated spins and 
small clusters towards the largest cluster, (a): Initially, a 
large square cluster is surrounded by 4 isolated "up" spins and 
3 smaller clusters (b): A shear step occurs and the top half of 
the simulation box is shifted by one lattice space to the right 

(c) : The cluster grows in the ±a; directions via Metropolis 
spin flips enhanced by the concave kinks created by the shear 

(d) : In the resulting configuration, the largest cluster has 
coalesced with one smaller cluster and has also become one 
lattice space closer to the small cluster at the bottom right. 



A second possible mechanism for the increase in / with 
7 is shear-induced coalescence between isolated up spins 
or small clusters and the growing nucleus. The shear al- 
gorithm is expected to drive together isolated up spins, 
causing an increased abundance of small clusters in the 
system, as illustrated in Figure [H These may then co- 
alesce with the largest cluster. Moreover, isolated spins 
and small clusters can also be driven towards the largest 
cluster by the shear algorithm. This is illustrated in Fig- 
ure [71 A shear step creates a kink in the largest cluster 
(FigEla — > b). This is a favourable site for growth, which 
tends to fill in the clefts created by the kink - with the 
result that the cluster grows preferentially in the x di- 
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rection (Fig[7]c). This growth reduces the gap between 
the largest cluster and surrounding clusters in the ±x 
directions (Fig[7]d). Alternatively, multiple shear steps 
occurring at the second row from the top in Figure [7h 
would shift the small cluster at the top left of the box 
towards the largest cluster, eventually resulting in coa- 
lescence, without the need for any growth of the largest 
cluster. Figures [6] and [7| demonstrate that our shear al- 
gorithm can promote growth of the largest cluster by 
coalescence, even though the system dynamics consists 
only of Metropolis spin flips and shear steps (diffusion of 
spins is not modelled). 



VI. ANALYSIS OF CLUSTER GROWTH 



To elucidate the role of enhanced cluster growth and 
coalescence in shear-enhanced nucleation, we first anal- 
ysed the ensemble of transition paths generated by the 
FFS simulations at shear rates 7 = 0.0, 7 = 0.06 and 
7 = 0.12. For each shear rate, we computed the con- 
tributions to the growth of the largest cluster of "single 
spin flip growth" (events in which the size of the largest 
cluster changes by ±1 spin) and of coalescence (events 
in which the size of the largest cluster increases by more 
than one spin). We also measured the contribution of 
cluster breakup events, in which the size of the largest 
cluster decreases by more than one spin (an "event" here 
refers to a single attempted spin flip in our Metropolis 
Monte Carlo scheme, or a shear step). The results are 
plotted as an average over the transition path ensemble, 
as a function of the committor Pb, in Figure [H The com- 
mittor function Pb{x) is the probability that a trajectory 
initiated from a conflguration x will reach the flnal state 
B before the initial state A - this provides a convenient 
measure for the progress of the transition. 

Figure [HK shows that, for all shear rates, the largest 
cluster increases in size as the transition progresses in a 
rather similar manner, although the shear causes a slight 
difference in cluster size: Nc{j = 0.06) < A^c(7 = 0.0) < 
^cii = 0.12). However, Figure [Hb, c, and d show that 
despite this apparent similarity, the contributions of sin- 
gle spin flips, coalescence and breakup events to clus- 
ter growth are all strongly affected by the shear. Both 
single spin flip growth and coalescence are enhanced for 
7 = 0.06 compared to the zero shear case, and strongly 
enhanced for 7 = 0.12. However, this is balanced by 
a strong increase in the negative contribution of clus- 
ter breakup events in the presence of shear. This sug- 
gests that both the mechanisms outlined above (shear- 
enhanced single spin flip growth and shear-enhanced co- 
alescence), as well as cluster breakup, are likely to play a 
signiflcant role in the nucleation mechanism in the pres- 
ence of shear. 
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Figure 8: (a): Largest cluster size Nc plotted as a function 
of the committor Pb, averaged over 25 transition paths, (b): 
Contribution of single spin flips {i.e. spin flips where ANc = 
±1) to Nc. (c) Contribution of coalescence events {i.e. spin 
flips where ANc > 1) to Nc. (d) Contribution of breakup 
events {i.e. spin flips where AA'c < — 1) to A^c. In all plots, 
circles represent results for 7 — 0.0, squares for 7 = 0.06 and 
diamonds for 7 = 0.12. Note the different scales on the Nc 
axis. 



VII. A MODIFIED SHEAR ALGORITHM 

We would like to quantify the contributions of shear- 
enhanced cluster growth and coalescence to the increase 
in nucleation rate with shear rate shown in Figure [H To 
do this, we devised a modifled shear algorithm in which 
shear-induced coalescence is largely eliminated. This is 
the same as the regular shear algorithm (see section [Ul, 
except that after each shear step, spins in the rows imme- 
diately above and below the "break point" are allowed to 
equilibrate with no shear for N^q Metropolis Monte Carlo 
steps. During this equilibration, all spins in the largest 
cluster and its immediate neighbours are held fixed, even 
if they lie in the rows mentioned above. The effect of 
this equilibration is that the largest cluster (and its im- 
mediate neighbours) is sheared as normal, while the sur- 
rounding "bath" of spins is maintained close to equilib- 
rium (there may still be small clusters in the bath which 
have broken off from the largest cluster) . It is not neces- 
sary to equilibrate all rows after a shear step. Our system 
has only nearest-neighbour interactions, and the effect of 
the shear is simply to change the juxtaposition of the two 
rows adjacent to the break point; all other rows remain 
unchanged during a shear step. We have used Neq = 5 
equilibration steps after each shear step (although in fact 
we find that N^q = 1 is sufficient). 

We expect this modified algorithm significantly to re- 
duce shear-induced cluster coalescence. The shear cannot 
form small clusters or drive surrounding spins and clus- 
ters towards the largest cluster, as illustrated in Figure 
[7l Coalescence events in which the shear directly merges 
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(c) Pb (d) Pb 

Figure 9: Contributions to the largest cluster size as in Figure 
El for the modified and regular shear algorithms for 7 = 0.06. 
Red squares: modified shear algorithm; black squares: regular 
shear algorithm (both for 7 = 0.06). Results are compared 
to those for the regular shear algorithm in the absence of 
shear, shown by the black circles, (a): Largest cluster size 
Nc- (b): Contribution of single spin fiips (i.e. spin fiips where 
ANc = ±1) to Nc. (c) Contribution of coalescence events 
(i.e. spin fiips where ANc > 1) to Nc- (d) Contribution of 
breakup events {i.e. spin fiips where ANc < — 1) to A'^c. Note 
the different scales on the Nc axis. 



clusters, as illustrated in Figure [H may still occur, how- 
ever, as may events where the largest cluster is broken 
up by the shear and subsequently re-coalesces. This type 
of modified shear algorithm might also be useful for in- 
vestigating the effects of coalescence for more complex 
systems with long-range interactions or off-lattice par- 
ticles: however, in these cases it would be necessary to 
equilibrate the whole bath after every shear step (or from 
time to time in the case of continuous shear) . If the nucle- 
ation mechanism proceeds via the coalescence of multiple 
large clusters (at large supersaturation) , this approach is 
unlikely to be useful. However, our simulations are at 
moderate supersaturation so that we believe nucleation 
occurs via the growth of a single large cluster. 

Figure [9] shows the contributions to cluster growth of 
single spin fiips, coalescence and breakup (as in Figure[8l), 
for the modified shear algorithm, compared to the regu- 
lar algorithm, for 7 = 0.06 (moderate shear). Transition 
paths were extracted from a FFS calculations: for the 
modified algorithm, we used 32 interfaces and Xb = 805, 
as the algorithm is computationally expensive. We veri- 
fied that this change of FFS parameters has no effect on 
the computed rate constant, since by A = 805 the paths 
are completely committed to nucleation. 

As expected, single spin fiip growth is enhanced by 
shear for the modified algorithm in the same way as for 
the regular algorithm (in Figure[9b, the data for the mod- 
ified algorithm with 7 = 0.06 overly those for the regular 
algorithm with 7 = 0.06). Turning to Figure [91:, we see 



that the coalescence contribution is significantly reduced 
for the modified algorithm (red squares) compared to the 
regular algorithm (black squares). However, some shear 
enhancement of coalescence still remains in the modi- 
fied algorithm, since the results for the modified algo- 
rithm with 7 = 0.06 (red squares) do not coincide with 
those for the regular algorithm in the absence of shear 
(black circles). We speculate that these remaining coa- 
lescence events involve clusters breaking up and then re- 
coalescing (these events are not suppressed by our mod- 
ified algorithm). Another possibility may be that the 
shear changes the largest cluster in such a way that it 
becomes more prone to attaching other clusters, even in 
an equihbrium bath. Figure |9li shows the contribution 
of cluster breakup. Cluster breakup is also partially re- 
duced in the modified algorithm compared to the regular 
algorithm. We expect the largest cluster to be stretched 
by the shear in the same way for both algorithms - so 
we might expect the contribution of cluster breakup to 
be the same for both algorithms. However, some coales- 
cence events result immediately in breakup, if the incom- 
ing small cluster is not well attached to the largest clus- 
ter. Suppressing these coalescence events in our modified 
algorithm may therefore also decrease the rate of breakup 
events. 




Y 

Figure 10: / versus 7 for h = O.OSfcsT, for "regular" shear 
(circles) and the modified shear algorithm (squares). The 
"regular" shear results are the same as in Figure [H 



Figure [TOl shows the nucleation rate / as a function of 
the shear rate 7 for the modified algorithm, as well as the 
original results from Figure [H The modified algorithm 
gives the same qualitative behaviour for 7(7), but the 
nucleation rate is reduced by a factor which increases 
with the shear rate. These results suggest that shear- 
induced interactions between the growing cluster and the 
surrounding sheared "bath" of spins play a significant role 
in enhancing nucleation for shear rates 7 > 0.03. 
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VIII. ANALYSIS OF THE TRANSITION STATE 
ENSEMBLE 

To further test the role of cluster coalescence in the 
nucleation mechanism, we have analysed the committor 
function Pb{x). This is the probability that a trajec- 
tory initiated from a configuration x will arrive in the 
B state before the A state. Surfaces in configurational 
space on which Pb takes a fixed value are known as the 
isocommittor surfaces. In particular, the transition state 
surface has Pb — 0.5. It is important to be clear about 
which configurations x we use to define the isocommit- 
tor surfaces. Because our system is driven by shear and 
its dynamics do not obey detailed balance, the path en- 
sembles for the forward and reverse transitions need not 
lie in the same region of configurational space, so that 
the transition state surfaces for the forward and reverse 
transitions need not be the same for a nonequilibrium 
problem Even for cases where detailed balance is 
obeyed, it is important to define whether the isocom- 
mittor surface is computed using configurations from the 
path ensemble or using Boltzmann-weighted configura- 
tions sampled from the entire phase space. We will study 
here the transition state surface for the forward paths: 
i.e. we will analyse the collection of configurations in 
the forward path ensemble with Pb — 0.5. These con- 
figurations are members of the transition state ensemble 
(TSE) for the forward transition. We will carry out this 
analysis for three different shear rates 7 = 0.0, 7 = 0.06 
and 7 = 0.12. The TSE configurations were extracted 
from the transition paths as described in Appendix [Bl 

Committor analysis can be used to test whether a cho- 
sen order parameter fj, is important in the transition 
mechanism, by computing the correlation between the 
value of ^ and the committor value, for configurations 
in the transition paths. If fi is found to be strongly cor- 
related with the committor, then it is likely that this 
order parameter captures (at least some of) the essential 
physics underlying the transition. For this nucleation 
problem, we know that the largest cluster size Nc is an 
order parameter that couples strongly to the committor - 
large clusters have a greater probability of continuing to 
grow than small clusters. In Classical Nucleation The- 
ory, it is assumed that Nc is the only important order 
parameter. However, in the presence of shear, other or- 
der parameters, which couple to the shear, must also be 
important. Since we have postulated that cluster coales- 
cence plays an important role in the nucleation mecha- 
nism in the presence of shear, we seek an order parameter 
/i that is coupled to coalescence. 

Figure fTD illustrates one such order parameter: the lo- 
cal density of up spins surrounding the largest cluster. 
The coalescence mechanism illustrated in Figure [7| is ex- 
pected to depend on the density of up spins in the regions 
to the right and left of the cluster, but not in the regions 
directly above and below. We therefore define local den- 
sities Px and Py as the density of up spins surrounding 
the largest cluster, located a distance < ^ or |y| < ^ 
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Figure 11: Order parameters for local density of "up" spins: 
px is the density of "up" spins, not in the largest cluster, in 
the shaded region in panel (a), while py is the density of "up" 
spins, not in the largest cluster, in the shaded region in panel 
(b). 



from the edge of the largest cluster respectively, where ^ 
is a cutoff distance. These regions are shown schemati- 
cally in Figure [TTk and b. If shear- induced coalescence 
is indeed important in the nucleation mechanism, we ex- 
pect Px but not Py to be correlated with the committor 
in the presence of shear. Neither px nor py is expect to 
correlate with the committor in the absence of shear. 

To test whether px and py are correlated with the 
committor, we extract TSE configurations (which have 
Pb — 0.5) from our transition paths, and look for nega- 
tive correlation between px or py and the largest cluster 
size Nc for these configurations. We know that Nc is an 
important order parameter for this transition. If another 
order parameter p is also important then we would ex- 
pect both Nc and p to determine the committor value. 
Specifically, TSE configurations can achieve Pb = 0.5 by 
having a large value of A^c but only a small value of p, 
or alternatively a small value of Nc may be compensated 
by a large value of p. This impHes that if p is positively 
correlated with the committor, then for the TSE config- 
urations, p should be negatively correlated with Nc- 

Figure[T2ltests whether this negative correlation exists. 
The main panels show scatter plots of px versus Nc for 
configurations in the TSE (each configuration is repre- 
sented by one dot), for ^ = 5 lattice sites. In the absence 
of shear (7 = 0.0), there is no significant correlation be- 
tween Px and Nc, but significant negative correlation is 
observed for 7 = 0.06 and 7 = 0.12. The correlation coef- 
ficient ("r-value") is larger for 7 — 0.12 than for 7 — 0.06, 
suggesting that the importance of coalescence increases 
with shear rate. The insets in Figure [12] show the same 
scatter plots for py. No significant correlation is present 
in the absence of shear (a) or when 7 = 0.06 (b). This 
supports our hypothesis that local density in the x but 
not the y direction is important in the presence of shear. 
For the highest shear rate, 7 = 0.12, we do observe sig- 
nificant negative correlation between Py and Nc, but the 
magnitude of the correlation coefficient is much less for py 
than for Px. Taken together, these results suggest that 
the shear-induced coalescence mechanism illustrated in 
Figures [6] and [7| does capture some of the essential physics 
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Figure 12: Scatter plots of px (main panel) or py (inset) versus 
size of largest cluster Nc for TSE configurations (defined as 
configurations in transition paths with Pb ~ 0.5), for (a): 
7 — 0.0, (b): 7 = 0.06 and (c): 7 = 0.12. Order parameters 
px and py were defined as above with ^ = 5 lattice sites. 95% 
confidence intervals for the correlation coefficient ("r value") 
between p and are given in each panel [50l |. 



of the nucleation mechanism. 

Since we beheve that shear also enhances "single spin 
flip" cluster growth, we would like to repeat the above 
analysis for an order parameter that measures the ten- 
dency of the cluster to grow. However, this has proved 
difficult. Order parameters based on measuring the num- 
ber of kinks in the perimeter of the largest cluster tend to 
covary with the size of the largest cluster, so that scat- 
ter plots of versus Nc are not an objective measure. 
These order parameters are also affected by the shape 
of the cluster, which is affected by the shear. Attempts 
to measure directly the propensity of the cluster to grow 
by single spin flips failed to produce signiflcant negative 
correlation, possibly because this propensity fluctuates 
greatly during cluster growth. 



IX. DISCUSSION 

In this paper, we have used Forward Flux Sampling 
(FFS) to calculate rate constants and transition paths 
for homogeneous nucleation in a sheared two dimensional 
Ising model. Our results show a striking nonmonotonic 
dependence of the nucleation rate / on the shear rate 7. 
We have investigated the physical mechanisms underly- 
ing this behaviour by analysing transition paths and the 
transition state ensemble, as well as by comparison to 
several modifled shear algorithms. Our conclusions are 
that the observed decrease in 1(7) for large 7 is due to 
shear-mediated breakup of the growing cluster, while the 
increase in 1(7) for small 7 is due to shear-induced cluster 
coalescence as well as to shear-enhanced "single spin flip" 



growth of the largest cluster. The contributions of shear- 
enhanced cluster coalescence and single spin flip growth 
appear to be of the same order of magnitude. 

Our analysis has been strongly influenced by work on 
(quasi-)equilibrium systems, in which the goal is often 
to identify the "reaction coordinate". This is the collec- 
tive coordinate which most strongly correlates with the 
committor for the transition. If the true reaction coordi- 
nate could be identified, it is believed that the transition 
could be coarse-grained into a one-dimensional diffusion 
problem, over a free energy barrier in the reaction co- 
ordinate, in the manner of Classical Nucleation Theory. 
Several methods have been proposed for identifying the 
optimal reaction coordinate [34, 51, 52|- We were mo- 
tivated by this approach in our attempts to find order 
parameters that correlate negatively with Nc in Section 
IVl However, it is not clear whether the reaction coordi- 
nate is such a useful concept for driven systems as it is 
for equilibrium systems. As discussed in section lVIIH the 
isocommittor surfaces for the forward and backward tran- 
sitions are not necessarily the same for driven systems, so 
the reaction coordinates may be different for the forward 
and backward processes. Moreover, the reaction coor- 
dinate is usually assumed to be a function only of the 
spatial coordinates of the system - whereas in systems 
that are not in local equilibrium, dynamical coordinates 
(momenta) may also be important. Even if a collective 
coordinate which correlated precisely with the committor 
could be found, it is not clear whether the complex dy- 
namics of the driven system could then be coarse-grained 
into a one-dimensional model. Our preliminary investi- 
gations suggest that it may be possible to construct a 
one-dimensional "CNT-like" model for Ising nucleation 
under shear ^along the lines of a toy model proposed by 
Cerda et a/ ji,[H3|. 

The two dimensional Ising model studied in this pa- 
per is an idealised system, and the algorithm used here 
to apply shear is far from realistic. In particular, our 
dynamics includes no spin transport and our shear algo- 
rithm imposes a Hnear "velocity gradient" on the system. 
In experimental systems, mass transport, hydrodynamic 
effects, etc, undoubtedly play a role. Effects of mass 
transport could be included using an Ising model with 
Kawasaki dynamics [45| instead of MetropoHs spin fiips 
as used here. This would require some reservoir of up 
spins to be provided. In the absence of driving, the choice 
of dynamical update rule has very significant effects on 
system behaviour [53|: a comparison between Metropolis 
and Kawasaki dynamics for nucleation under shear would 
be an interesting topic for further work. Despite its sim- 
plicity, the Ising model has made important contributions 
to the understanding of nucleation phenomena in equi- 
librium systems, and is likely to yield important insights 
also for nonequilibrium systems. Moreover, the mecha- 
nisms identified here - shear- induced cluster breakup, en- 
hanced cluster growth and cluster coalescence - are Hkely 
to play a significant role in experimental nucleation under 
shear. The discreteness of our model may be a cause for 
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concern: for example, the "kink" mechanism for shear- 
enhanced single spin flip growth illustrated in Figured]) 
is a discrete phenomenon. However, our results show 
that cluster coalescence is of approximately equal impor- 
tance, and, in any case, preliminary results indicate that 
the observed qualitative behaviour is not dependent on 
the size of the largest cluster jH3| ■ We are encouraged by 
recent two dimensional Brownian Dynamics simulations 
of charge stabilised and attractive colloids, in which simi- 
lar nonmonotonic behaviour of the crystal nucleation rate 
was observed [§| , although the underlying physical mech- 
anisms were postulated to be somewhat different. Future 
work will investigate the role of transport processes for 
nucleation under shear in simple model systems and ex- 
tend our work to more complex models. 
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Appendix A: THE SHEAR ALGORITHM 

Here, we describe the algorithm used to apply shear 
to the Ising system. The algorithm is similar to that of 
Girillo et al [47| but rows are shifted by only one lattice 
site. The algorithm required the parameters Ms (number 
of attempted row shifts per row per MG cycle) and Ps 
(probability of a row shift happening in each attempt) 
to be assigned, for a given shear rate 7. This is done as 
follows: if 7 < 1, Ms = 1 and Ps=i. If 7 > 1, then Ms 
and Ps are allocated such that Ms is the smallest integer 
value such that MsPs — 7 and < Ps < 1- Although this 
is the strategy used to assign Mg and Ps in this paper, 
we have verified that our results do not depend on the 
choice of Mg and Ps for a given 7. 

We then simulate the system using the usual Metropo- 
lis Monte Garlo algorithm, except that after each MG 
cycle, we carry out Ms x L attempts to shear the system 
by one lattice site. Each attempt consists of: 

(1) choose a random number < r < 1 

(2) If r < Ps, go on to (3) 

(3) choose a row jc at random (using a second random 
number) 

(4) for all spins with y co-ordinate j > jc, carry out the 
X co-ordinate transformation i — > i 1 - taking proper 



account of periodic boundary conditions 

(5) update the stored information on y-direction periodic 

boundary conditions, as detailed below 
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Figure 13: Schematic illustration of a row shift in the shear- 
ing algorithm. Columns are labelled by letters. The square 
simulation box is shown flanked by the bottom and top rows 
of the neighbouring boxes in the up and down direction in the 
periodic array (shown in grey). The small arrows indicate the 
neighbours of each lattice site in the top and bottom rows of 
the box in the up and down directions. A shear step results 
in a shift of all rows in the top half of the box by one lattice 
site. This results in a change of identity of the neighbours 
of the top and bottom rows of the lattice, as shown by the 
slanting arrows. 

In the absence of shear, we apply periodic boundary 
conditions (PBC) in the x and y directions. In the pres- 
ence of shear, PBG apply in exactly the same way in 
the X direction, but care is required in the y direction. 
As illustrated in Figure [131 a shear step results in a 
displacement of the identity of the "up" neighbours of 
the top row and the "down" neighbours of the bottom 
row. Therefore, in addition to storing the spin vari- 
ables a at each lattice site, for every configuration, we 
also store two integers, gup and gdown- The identity of 
the "up" neighbour of the spin at position [j; L] (where 
rows are numbered 1 ^ L) is [i -I- gup — nL; 1], while 
the identity of the "down" neighbour of the spin at po- 
sition [i; 1] is [i + gdown — mL;L]. The integers n and 
TO are chosen to ensure that 1 < i + gup — nL < L and 
^ < i + gdown — niL < L. Each time a shear step is car- 
ried out, gup and gdown are updated by gup ^ gup - ^ 
and gdown gdown + 1- It is essential when carrying out 
the FFS simulations and when reconstructing the tran- 
sition paths from the FFS sampling that gup and gdown 
are stored for each configuration. 

Appendix B: SAMPLING THE TRANSITION 
STATE ENSEMBLE 

To obtain configurations from the transition path en- 
semble with Pb = 0.5 (TSE configurations), we first re- 
generate the transition paths from the FFS simulation, as 
described in section lUl Having regenerated a transition 
path, we then fire 16 trial runs from every 10th configu- 
ration along this path, and monitor the number of trials 
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iVi^^ that reach B before A. If 7 < ivj^^ < 9, we then 
fire a further 100 trial runs from this configuration, and 
monitor the number trials Ns that reach B rather than 

(2) 

A. If 40 < Ns < 50, we consider this configuration a 
member of the Pb — 0.5 ensemble. This procedure pro- 
duces a collection of configurations with a range of Pb 



values around the desired value of 0.5. The parameters 
of the method can be adjusted to balance computational 
expense with accuracy. This approach is slightly differ- 
ent from that of Pan and Chandler [s^ , who include both 

Ns^^ and A^i^"* in the final computed value of Pb- 
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